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SECTION t 


INTRODUCTION 

This report summarizes the results of the work performed under NASA 
contract NASI -18098 from September 1985 through August 1986. The report is 
divided into five sections. Section 2 presents a new class of closed-form 
solutions for finite-time linear-quadratic optimal control problems, which 
is shown to be computationally more efficient than previously known closed- 
form solutions. Section 3 utilizes the closed-form solutions of Section 2 
for the feedback gains in the free-final-time perturbation feedback 
problem, where the initial conditions and terminal constraints may be 
assigned off-nominal values. Section 4 presents a control scheme for 
general nonlinear three-axis slewing maneuvers of flexible spacecraft. 
Under this control scheme, an open-loop rigid body nominal solution is 
applied to the spacecraft while a perturbation feedback controller reduces 
the elastic response and causes the system to closely follow the nominal 
rigid body trajectory. A modified Kalman filter is implemented for 
estimating the states of the system. Section 5 presents a summary and 
conclusions for this report. Reference 9 documents the detailed 
derivations of results presented in this report. 


I 
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SECTION 2 


CLOSED-FORM SOLUTIONS FOR FINITE-TIME LINEAR-QUADRATIC 
OPTIMAL CONTROL PROBLEMS 


2.1 Introduction 


During the design and analysis phases of optimal control synthesis, 
state and control trajectories are often computed to help the control 
engineer evaluate the control design. The most straightforward and most 
widely practiced method of computing the state and control trajectories is 
by numerically integrating the governing differential equations. This may 
be costly for flexible space structures which may have many elastic degrees 
of freedom. The reason for the high cost is two-fold: first, the large 
number of elastic degrees of freedom requires a large number of states to 
be integrated; and second, since the highest frequency of the system to be 
simulated increases as the number of elastic modes is increased, the 
integration step-size must be decreased correspondingly. Thus, the 
computational cost of the simulation increases rapidly as more elastic 
degrees of freedom are included. 

This section presents a new class of closed-form solutions for 
finite-time linear-quadratic optimal control problems when the plant is 
time-invariant. With a closed-form solution, one can compute the response 
of the entire system at any point in time. Thus, the engineer can compute 
the system response at any desired interval, independent of system 
frequency. (Of course, if one needs time-history plots with good 
resolution, the time interval does depend on system frequency.) In 
addition, numerical roundoff errors are greatly reduced, because the number 
of floating point operations 1 (flops) required for computing closed-form 


1 A floating point operation is more or less the amount of work needed to do 
a floating point add, a floating point multiply, and a little subscripting. 
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j 

solutions is much less than for numerical integration. Sensitivity 
partials may also be computed easily when closed-form solutions are 
available [35]. 

Other forms of closed-form solutions exist. For example, classical 
state and co-state solutions can be obtained from matrix exponential 
solutions where the system Hamiltonian matrix is used. However, such 
solutions either suffer from numerical instability or require too much 
computation. With a numerically stable and efficient closed-form solution 
and with the rapid development of parallel processing technology, one may 
envision the day when feedback gains may be computed on-orbit in real time. 
Whether the solutions proposed in this chapter can fulfill such a goal is 
of great interest, and is a subject for further research. 

Although there are many different finite-time optimal control 
problems, their solutions can be written as differential equations of only 
a few basic forms. The closed-form solutions of these basic differential 
equations are presented in Section 2.2, and example applications are 
provided in Section 2.3 for illustration. 

2.2 Closed-Form Solutions of Basic Differential Equations 

The solution of finite-time linear-quadratic optimal control 
problems involves the solution of differential equations which may be 
classified into five basic types: 

i 

j Type 1 


P(t) - - P( t) A - A T P(t) + P(t)EP(t) - Q , 

Type 2 


( 2 . 2 . 1 ) 


^(t) - - [A - EPUn^U) + F 1 (t) , 


( 2 . 2 . 2 ) 
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Type 3 


X 2 (t) = [A - EP(t)]X 2 (t) + F 2 (t) , 


(2.2.3) 


Type 4 


Y^t) 


X 1 1 a (t)EX 1b (t) , and 


(2.2.4) 


Type 5 


V t} = X 2a (t)Z 1 (t)EZ _1 (t)x 2 b (t) 


(2.2.5) 


Type 1 represents the well-known differential matrix Riccati 
equation with constant coefficients. Its solution, P(t), couples into the 
differential equations of Type 2 and Type 3- The functions X 1a (t) and 
X 1b (t) are solutions of differential equations of Type 2, and the functions 
X 2a (t) and X 2b (t) are solutions of differential equations of Type 3. In 
the above equations, A, E, and Q are (n * n) constant matrices, and the 
variables have the following dimensions: 


P(t) 

(n 

x n) 

x,(t) 

(n 

x p) 

X 2 (t) 

(n 

x q) 

Y 1 ( t) 

(r 

x s ) 

X 1a (t) 

(n 

x r) 

x 1b ( t) 

(n 

x s) 

Y 2 (t) 

(1 

x m) 

x 2 a (t) 

(n 

X 1) 

x 2 b (t) 

(n 

x m) 

The matrices E, Q, 

and P(t) 

must 


is a term representing the forcing functions for the X^(t) differential 
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equations, and accordingly has dimension (n x p) for i = 1 and (n x q) for 
i = 2. For scalar control problems, the five types of differential 

equations become scalar differential equations. However, for multivariable 
control problems, P(t) is always a square matrix; ( t ) (i-1,2) may 

represent either a matrix or a vector; and Y i ( t ) (i-1,2) may be either a 
matrix, vector, or scalar. 

2.2.1 Solution for Type 1 Differential Equations 

The Type 1 differential equations are defined by 

P(t) = - P(t)A - A T P (t) + P(t)EP(t) - Q . (2.2.6) 

The solution of the above differential matrix Riccati equation is well- 
known and, in fact, can be expressed in several different forms. One of 
the most useful forms of the solution is due to Potter [5,17,22,26-28], and 
is given by the sum of the steady-state solution (i.e. the solution to an 
algebraic Riccati equation) [1,5,17,22,25], and a transient term: 

P(t) = P ss + Z-V) , (2.2.7) 


where 


0 = P A + A P -PEP + Q . 
ss ss ss ss 


In order to obtain the differential equation for Z(t), one needs 
the following expression for the derivative of a matrix inverse: 



Substituting 


Z(t) 


V)] = - Z -1 (t)Z(t)Z -1 (t). (2.2.8) 

(2.2.7) into (2.2.6), and making use of (2.2.8), one obtains 
= AZ(t) + Z(t)A T - E , (2.2.9) 
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where 


A = 


A 


EP 


ss 


is the closed-loop system dynamics matrix. The solution for Z(t) can be 
cast in either of the following two forms [11,17,31]: 


A(t-t ) 

Z(t) = Z =c + e [Z(t )-Z ]e 

ss o ss 


A T (t-t Q ) 


, Z(t o ) 


[P(t )-P ] 

o ss 


-1 


or 


-A(t f -t) -A T (t f -t) . 

Z(t) = Z ss + e [Z(t f )-Z gs ]e , Z(t f ) = C p (t f )-P ss ] ; 


( 2 . 2 . 10 ) 


where Z ss satisfies the algebraic Lyapunov equation [2,13,17,24]: 


0 



+ Z 



- E . 


It can be shown that Z(t) and Z ^(t) exist for well-posed optimal 
control problems. 


As shown in Sections 2.2.2 - 2.2.5, the symmetric matrix Z(t) plays 
a central role in the solution of differential equations of Types 2, 3, 4, 
and 5. 


2.2.2 Solution for Type 2 Differential Equations 

The Type 2 differential equations are characterized by differential 
equations with time-varying coefficient matrices, where the coefficient 
matrices are functionally dependent upon the Type 1 equations. The general 
form for the Type 2 equations is given by 
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( 2 . 2.11 ) 


X^t) = - [A - EP(t) ] T X 1 (t) + F 1 (t) . 

On assuming a solution of the form 
X^t) = Z~ 1 (t)W(t) , 

where W(t) is unknown, it can be shown that the solution of 
given by [9] 


X^t) = 4> 1 (t,t o )X 1 (t o ) + 


t <t 1 (t,T)F 1 (t)dT , 


where 


-i X(t-t 0 ) 

$ 1 (t,t o ) = Z (t)e Z(t 0 ) 

is the state transition matrix for the homogeneous part of (2. 
integral term in (2.2.13) is easily obtained when F^t) can be 
terms of exponential matrices. 

2.2.3 Solution for Type 3 Differential Equations 

Like the Type 2 differential equations, the Type 3 
equations are characterized by differential equations with 
coefficient matrices, where the coefficient matrices are 
dependent on the Type 1 equations. The general form for 
equations is given by 

X 2 (t) - [A - EP(t)]X 2 (t) + F 2 (t) . 


( 2 . 2 . 12 ) 
(2.2.11) is 


(2.2.13) 


2.11). The 
expressed in 


differential 
time- varying 
functionally 
the Type 3 

(2.2.14) 


-T- 


On assuming a solution of the form 


x 2 (t) = Z(t)W(t) , 


(2.2.15) 


where W(t) is unknown, it can be shown that the solution of (2.2.14) is 
given by 


x 2 (t) 


* 9 (t,t )X_(t ) 

2 0 2 O 




4> 2 (t,T)F 2 (T)dT , 


( 2 . 2 . 16 ) 


where 


-A T (t-t ) 

* 9 (t,t ) = Z(t)e Z (t ) 

2 0 O 

is the state transition matrix for the homogeneous part of (2.2.14). The 
integral term in (2.2.16) can be easily evaluated when F 2 (t) is expressed 
in terms of products of Z(t) and exponential matrices. 

2.2.4 Solution for Type 4 Differential Equations 

The Type 4 differential equations are characterized by products of 
Type 2 solutions. The general form for the Type 4 equation is given by 

V t} “ X la (t)EX 1b (t) ’ (2.2.17) 

where the equations for X-j a (t) and X^(t) are, in general, inhomogeneous: 

X 1ft (t> - - [A - EP(t)] T X 1a (t) + F la (t) , (2.2.18) 

and 

X lb (t) = - [A - EP(t)] T X 1b (t) + F lb (t) . (2.2.19) 

From Reference 9, the solution for Y i ( t ) is shown to be 
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V‘> - x i T a< t ) z(t)5[ i 6 ' t > * W ' X 1a (t o )Z<t o ,X lb (t o ) 


- ft X t T a (,)Z(,)F lb (,)dT - ft F 1 T a (T)Z(T)X lb (T)<J, 


( 2 . 2 . 20 ) 

The integral terms in (2.2.20) are easily computed when F 1a (t) and F lb (t) 
are functions of exponential matrices. 

2.2.5 Solutions for Type 5 Differential Equations 


The Type 5 differential equations are characterized by products of 
Type 3 solutions and Z~^(t). The general form for the Type 5 equations is 
given by 


Y 2 (t) = X 2 T a (t)Z _1 (t)EZ 1 (t)X 2b (t) , 


( 2 . 2.21 ) 


where the differential equations for X 2a (t) and X 2b (t) are, in general, 
inhomogeneous: 

X 2a (t) = [A - EP(t)]X 2a (t) + F 2a (t) , (2.2.22) 

and 

X 2b (t) = [A - EP(t)]X 2b (t) + F 2b (t) . (2.2.23) 


From Reference 9, the solution for Y 2 (t) is shown to be 

Y (t) = - X* (t)Z -1 (t)X-.(t) + Y.(t ) + X* (t )Z _1 (t )X 

2 2a 2b 2 o 2a o o 


2b 


(t o> 


t 4 a ( T)Z- , (T ) F 2b (,)dt * 
0 


t f , 2a (T)Z " 1(l)X 2b (T)dT 
o 
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( 2 . 2 . 211 ) 


The integral terms in (2.2.2k) are easily computed when F 2 a (t) and F 2 b(t) 
can be expressed as products of Z(t) and exponential matrices. 

Throughout the developments of this section, one can observe the 
close relationship between equations of Type 2 and Type 3, and also between 
equations of Type 4 and Type 5. Indeed, this follows because equations of 
Type 2 and Type 3 are formal adjoints of one another. 

Tables 2-1 and 2-2 provide a summary of the five basic differential 
equations and their solutions. 

2.3 Example Applications of Closed-Form Solutions 


Reference 9 presents solutions of three finite-time linear- 

quadratic optimal control problems using the closed-form solutions of 
Section 2.2. For comparison, alternative closed-form solutions based on 
the state transition matrix of the state-costate system are presented. The 
comparison of the amount of computational work required for each type of 

solution clearly demonstrates that the new class of solutions is more 

efficient. 

In particular, using Potter's solution of Section 2.2.1, the 
propagation of the Riccati matrix is written as 

P(t + At) = P oa + Z"\t + At) , (2.3.1) 

SS 

where Z(t + At) is computed via 

Z(t + At) = C + e AAt Z(t)e A At , C = Z - e AAt Z e A At . (2.3.2) 

SS So 

Equation (2.3.2) requires roughly 3n^/2 flops for the propagation of the 
symmetric Z(t) matrix. In computing the Riccati solution of (2.3.1), the 
symmetric definite matrix inversion requires n^/2 flops. Thus, a total of 
2n^ flops are required to propagate the Riccati solution over one time- 
step . 
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Table 2-1. Summary of Basic Differential Equations and Their Solutions 



■ 11 ' 















An alternative way of computing the matrix Riccati solution is the 
well-known Kalman-Englar method [18], where P(t) is propagated at intervals 
of At by 


P( t+At) = [0 21 (t+At,t) + 6 22 (t+At,t)P(t)][0 11 (t+At,t) + 0 12 (t+At,t)P(t)] 

(2.3.3) 

where ©^(t+At.t) are partitions of the transition matrix for the state- 
costate system; that is, 


0 ( t + At, t) = 0( At) - e fiAt (2.3.4) 

and 

fi = 

For linear time- invariant systems, 0 is only a function of At, and hence 
need only be computed once. 

The number of operations required for the propagation of P(t) via 
(2.3.3) is n^ flops for each of the 0 2 2 P and e i2 p Products, n^/3 flops for 
the L-U decomposition (with partial pivoting) of the [0^ + 0 12 P] term, 
n^/2 for forward elimination, and n3/6 for back-substitution, where 
symmetry of P(t) is taken into account. The total number of operations 
adds up to about 3n^ flops. Thus, the Kalman-Englar method requires about 
50% more operations than Potter’s method for propagating P(t) at intervals 
of At. Moreover, numerical difficulties arise in the Kalman-Englar method 
when At is chosen too large, causing the term to be inverted to be nearly 
singular [21]. For the propagation based on Potter’s solution, such 
difficulties do not occur. Another solution approach is the negative 
exponential solution derived by Vaughan [ 38 ], which produces a numerically 
stable algorithm. However, since this method involves complex eigenvectors 
of the Hamiltonian matrix, the use of complex arithmetic causes the 
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operation count to be many times higher than that for the Kalman-Englar 
method. 


Table 2-3 summarizes the solutions for state and control 
trajectories of the three example problems presented in Reference 9, which 
include the optimal linear regulator, the controller with terminal 
constraints, and the tracking/disturbance accommodating controller. A 
remarkable fact is that despite the differences between the three control 
problems, the state and control trajectories may be cast into the same 
general form, namely, 

x ( t ) = Z a(t) + b(t) , (2.3.5) 

ss 

and 

u(t) - - D 1 a(t) - D 2 b(t) , (2.3.6) 


where 


1 B T 


[p z 

L S3 SS 


1] 


and 


•1 T 
BP 


ss 


The definitions of a(t) and b(t) , however, are slightly different for each 
problem. 


2.H Subspace Reduction for the Hamiltonian Matrix 

It can be shown that the new class of solutions, which involves the 
variables P ss , Z QQ , e^ fc , and e~^ Tt , are related to the closed-form 
solutions involving partitions of the state transition matrix by means of 
reducing subspace transformations of the Hamiltonian matrix: 


Table 2-3. Summary of Example Applications of Closed-Form Solutions 


























Table 2-3. Summary of Example Applications of Closed-Form Solutions (Continued) 



tion to x(t) 
and u(t) 
(flops) 


















(2.4.1 ) 


A 

- E 1 


i 

Z ss 


■ A 

0 


' I+Z ss P ss 

■ Z ss" 

. " Q 

1 

> 

-3 

1 


- P ss 

P Z +1 

SS SS -J 


0 

-; T . 


- - p ss 

I 


In (2.4.1), the left hand side represents the Hamiltonian matrix, ft, of 
(2.3.4). Exponentiation of (2.4.1) leads to the following block- 
diagonalizing transformation for the state transition matrix: 


A(t-t 0 ) 


r o.. 0. 


' I Z 


■ e 0 o 

11 12 


SS 


-T 





-A X (t-t ) 

0 0 


P P Z + I 


n 0 

U 21 22 J 


L SS SS SS J 


L 0 e J 


X 


I + Z P 

S3 S3 


-P 

SS 



(2.4.2) 


where 0^ = 0jj(t,t Q ) are partitions of the state transition matrix 


ft(t-t Q ) 

0(t,t Q ) = e 


(2.4.3) 


-IT- 


SECTION 3 


SPACECRAFT SLEWING MANEUVERS USING A CLOSED-FORM SOLUTION 
FOR THE NEIGHBORING EXTREMAL PATH PROBLEM 


3.1 Introduction 


The optimal control problem in this section is specified by 
defining a performance index which consists of a penalty on elapsed time, 
a quadratic penalty on the terminal states and controls, and an integral 
of quadratic penalties on the states, controls, and control rates. The 
final time is free, and specified terminal constraints produce a terminal 
manifold which also may be a function of the final time. Assuming that 
the nominal control and state trajectories are known, one seeks the 
perturbation feedback gains which cause the system to follow a neighboring 
extremal path when subjected to small perturbations in the initial 
conditions and terminal constraints. Necessary conditions for the 
perturbed system are stated, and the solution for the nominal trajectory 
is shown. Solutions for the perturbation feedback gains are developed 
based on the results of Section 2. Perfect plant knowledge and perfect 
state estimation is assumed. A time-to-go indexing scheme is used for 
applying the feedback gains so that the controller does not run out of 
feedback gains if the actual final time is longer than the nominal final 
time. Slight numerical modifications are presented for overcoming the 
numerical sensitivities of this type of controller. Two retargeting 
example maneuvers are shown, involving a spacecraft model consisting of a 
rigid body with four flexible appendages. An extension is proposed for 
using the closed-form solutions in control problems involving nonlinear 
systems by linearizing the nonlinear plant equations about the nominal 
trajectory. 
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3.2 


Statement of the Control Problem 


Let us assume that we have obtained the p-dimensional nominal 
control vector u^(t), which minimizes the quadratic performance index 


J 


Vr * ?fVf 




x^W x + u*W u]dt , 
xx uu 


( 3 . 2 . 1 ) 


subject to 

x = Ax+Bu, x(t) = x given , (3.2.2) 

o o 

N 

<|>Lx(t ),t ] = Mx(t f ) - ^ (t f ) = 0 , and t f unspecified . 

(3.2.3) 

In the above equations, x is the n-dimensional state vector, A and 

B are the time-invariant state dynamics and control influence matrices, 

T 

is a q-dimensional vector of terminal constraints, > 0 and W xx = 

W xx > 0 are weighting matrices for the state, W uu = W^ u > 0 is a weighting 
matrix for the control, and W t > 0 is a weight for the final time. The 
following necessary conditions must be satisfied by the nominal optimal 
control and state trajectories [7]: 


& = Ax+Bu , x(t Q ) given , (3.2.4) 

X = - W xx x - A T A , A(t f ) = S f x(t f ) + M T v , (3.2.5) 

u = - W -1 B T A , (3.2.6) 

uu 

<|<[x(t f ) ,t f ] = 0 , (3-2.7) 
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where 


$ = W t fc f + 1 x f S f x f + vT * » 

f - w t - ^d ( V + [x r s f + v ' M]x f ’ 

A ( t ) is the n-dimensional costate vector, and v is the q-dimensional 
vector of Lagrange multipliers for the terminal constraints. 

Let us now consider small perturbations in the initial states 
5x(t 0 ), and in the terminal constraints dip . The perturbation problem is 
then to seek the the correction to the control, 6u(t), which causes the 
perturbed system to minimize the original performance index subject to the 
new initial conditions and new final constraints. Moreover, we seek a 
feedback form for the solution of 6u(t), which involves the perturbations 
in the state, <Sx(t), and the perturbations in the final conditions dip . 
The necessary conditions which must be satisfied by the perturbed system 
are given by the following equations [7]: 


6x = A6x + BSu , 6x(t Q ) given , 
6A = - W xx 6x - A T 6A , 


(3-2.9) 

( 3 . 2 . 10 ) 


6u = - wJb T SA , (3.2.11 ) 

•»<v • [75-11 8x f . [|f]^ dv * [f^r dt f 

3x t f t f t f 

= S f 6x + M T dv + [S f x f + A T S f x f + A T M T v + W xx x f ]dt f , (3.2.12) 

d<l> = [f£]| 6x f + [-^]| T dt f = M6x f + [Mx f - ^(t f )]dt f , 
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and 


« - [|f]| sx f * [f]| dv * [f]| dt f - 0 


= [xjs f + xJs f A + v T MA + xjw xx ]6x f + [Mx f - \|/J(t f )] dv 


[-v\(t f ) + (xjs f + x*S f A + v X MA + xJ.W xx )x f ]dt f 


(3.2.14) 


3.3 Solution for the Nominal Trajectory 


The solution for the nominal trajectory may be obtained by using 
the state transition matrix and an exponential form for J. For a given 
final time tf, the final states and costates can be written as 


and 


' x(t f ) ' 


G(t--t ) 
r f 01 


x(t o } I 

. . 

► = 

[e J 

< 

► 

A (t ) 

O' 


where 


A 

-W 


xx 




°<V t o ) 


xx 


AX 


xA 


AA J 


(3.3.1) 


is the state transition matrix for t * t^. By introducing the terminal 
constraint of (3-2.3) into (3.3.1), and using the boundary condition of 
(3.2.5) , one obtains 
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I 

s. 


0 


M 


f 

M 0 


' x(t f ) ' 


'T 

XX 

XA 


X( V I 

4 

V 

► = 

*Ax 

A 

^AA 

A 

4 

► 

A (t ) 

0 


► + i 


0 

0 


I J 

(3.3.2) 


Upon rearranging and placing the unknown variables on the left hand side, 
one obtains 


■ I 

0 

, 1 
xX 


' x(t f ) ' 


$ 

XX 

0 “ 


s f 

T 

M 

XX 

< 

V 

► =s 

4> 

Xx 

0 

* 

M 

0 

0 


. »<t 0 > . 


0 

I _ 



♦S'V 


^ ,(3.3.3) 


from which x(tf), v, and A(t 0 ) may be obtained via Gaussian elimination. 


The above solution assumes that tf is known. However, for a free- 
final-time problem, the optimality condition, fi = 0, of (3.2.8) must be 
satisfied as well. This conditon produces a local minimum for J. To 
obtain the global minimum, one can numerically compute the value of J over 
a reasonable range of t^ and find which value of t^ produces the lowest 
value of J. Efficient propagation algorithms for J are presented in 
Reference 9 for computing values of J for different final times. 


Using (3.2.1), (3.2.8), (3.3-3), and the propagation equations for 
J, one can numerically compute values of J and ft over a range of final 
times, and hence find the optimal t^. For the case where one wishes to 
adjust Wj. so that the optimal final time is at a desired value, one merely 
computes a using (3.2.8) and (3.3-3) at the desired final time and obtain 
the required value of to make U = 0. 

Having found the values of the optimal final time and initial 
costates, the nominal state and control trajectories are given by 
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and 


N,. v 
x (t) 


= [I 0]e 


G(t-t ) 
o 


x(t o } 

x(t o ) 


N,. . 
u (t) 


G(t-t ) 

W V[0 I]e 

uu 


x(t o ) 

A(t 0 ) 


( 3 . 3 . 4 ) 


( 3 . 3 . 5 ) 


3.4 Solution for the Feedback Gains 


We now seek the solution for 6u(t) in the form 


6u(t) = - K 1 (t)6x(t) - K 2 ( t )d^ , 


(3.4.1) 


where 


u(t) = u N (t) + 5u(t) , 


6x(t) = x N (t) - x(t) , 


and K 2 are the required feedback gains, and u(t) and x(t) are the 
perturbed controls and states. By manipulating the terminal conditions of 
(3.2.12) through (3.2.14) into the following form, which is assumed to be 
valid for t Q <, t £ t^, the costate perturbations are expressed in a 
feedback form which leads to (3.4.1): 


' 6A(t) ' 


‘ S(t) 

R(t) 

m(t) " 


' <5x(t) ' 

dv 

► = 

T 

R (t) 

S(t) 

n(t) 

1 

—dip ► 

dt f 


_ m (t) 

n (t) 

a(t) _ 


-dn = 0 

\ 


(3.4.2) 


- 23 - 


where 


S(t f ) = S - RQ 1 R T , R(t f ) = - RQ _1 , m(t f ) 

Q(t f ) = - Q _1 , n(t f ) = Q _1 n , S(t f ) = 



and 



Notice that Q is singular because Q = 0. 


--- 1 _ 

RQ n - m , 
n Q n - a , 


Treating dv, dtf, d\|>, and dft as constants, one can differentiate 
(3.4.2) using (3.2.9) through (3.2.11) and collect terms to obtain the 


differential 

equations for the 

coefficient matrices: 


• 

S 

= ■ 

A T~ A AS 

- SA - A S + SES - 

W xx * 

(3.4.3) 

• 

R 

* 

- [A - ES] T R , 


(3.4.4) 

• 

A 

m 

= 

- [A - ES] T m , 


(3.4.5) 

• 

/V 

Q 

= 

aT /a 

R ER , 


(3.4.6) 

• 

n 

rs 

R Em , 


(3.4.7) 

• 

a 

=s 

A T -. 

m Em , 


(3.4.8) 
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where 


E = 


-1 T 
BW B 
uu 


A /V 

The matrices S(t) and R(t) are used for the perturbation feedback, 
while the vectors ffi(t) and A(t) are used for the estimation of the change 

A 

in the final time. The matrix Q(t) and the scalar S(t) are not needed for 
solving the control problem. 

A 

The solution of the Type 1 differential Riccati equation for S(t) 
is shown in Section 2.2.1 to be 

S(t) = S + Z _1 (t) , (3.4.9) 

ss 

A 

where S gg satisfies the steady-state Riccati equation 


0 = - 


S A - 
ss 


T~ 

AS 


ss 


+ SES - W 


XX 


From Table 2-2, the solution for Z(t) can be shown to be 
A(t-t ) A T (t-t ) 

Z(t) = Z + e ° [Z(t ) - Z ]e ° , (3.4.10) 

SS o ss 

— A 

where A = A - ES gg and Z gg satisfies the steady-state Lyapunov equation 

for Z(t) defined following (2.2.10). Since R(t) and in(t) satisfy differ- 
ential equations of Type 2, their solutions are given by (see Tables 2-1 
and 2-2): 


A( t-t ) 

R(t) = Z (t)e Z(t o )R(t Q ) , 


(3.4.11 ) 


and 


_ A(t-t ) 

fn(t) = Z (t)e Z(t )m(t ) . 

o o 


(3.4.12) 
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Since n(t) satisfies a differential equation of Type 4, its solution is 
obtained from Table 2-1 as 

n(t) = R T (t)Z(t)m(t) + n(t Q ) - R T ( t Q ) Z ( t Q ) fh( t Q ) . (3.4.13) 

The initial conditions for S(t), R(t), m(t), and n(t) are required for 
propagating the feedback gains forward in time; these are listed in 
Reference 9. 

Using (3.2.11) and (3*4.2), one can write the feedback gains of 
(3.4.1 ) as 

K 1 ( t) - W~Vs(t) , ( 3 . 4 . 14 ) 

and 

K 2 (t) = “ W uu Bl * (t) * (3.4.15) 


where the solutions of S(t) and R(t) are given by (3.4.9) and (3.4.11), 
respectively. 

3.5 Time-To-Go Indexing Scheme 


On observing the time arguments in (3.4.1) and recalling the fact 
that the final time is free, it quickly becomes apparent that if the 
optimal final time of the perturbed system lies beyond the nominal final 
time, then the feedback gains are undefined for part of the time (t > t^) 
along the neighboring extremal path. One of the methods suggested for 
eliminating this problem is the use of time-to-go indexing [3,20,33] so 
that the time-to-go on the perturbed trajectory is the same as the time- 
to-go on the nominal trajectory (see Figure 3-1). Equation (3.4.1) is 
then re-written as 

6u(t) = - K 1 (t I )6x(t) - K 2 ( tjjdip , (3.5.1) 
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t = Current Time 


Indexed Time 


t 


t = Actual Final Time 
Gains are Computed 


Figure 3-1. Time-To-Go Indexing Scheme 
-21 


t * = Nominal Final Time 
t = Time-points at uih i ch 



where 



t ™ t j , 

N 

t is the current time, tj is the indexed time, tf is the nominal final 
time, and t f is the final time on the perturbed trajectory. 

To compute the indexed time, let us re-write the last row of 
(3.4.2) for the change in final time: 

dt f = m T (tj)6x(t) - n T (tj)di|j , (3.5.2) 


'"‘T -T 

where the time arguments for m (t) and n (t) have been changed to tj. 

A A A A 

Since the feedback gains S(t), R(t), m(t), and n(t) are efficiently 

propagated at fixed time intervals, let us define tj as the points in 
time at which the values of the feedback gains are available: 

tj = nAt, n = 1,2,3, ... (3*5.3) 

where At is the propagation time step. On assuming that the perturbed 
terminal manifold is given by 

Mx(t f ) - W , (3.5.4) 


the vector d\p of (3.5.2) is computed via 
dip = 'i'dttf ) "* ) • 


(3.5.5) 


Since t f of (3.5.5) depends on the value of dt f , (3.5.2) represents an 
implicit equation for dtf. 

To interpolate the gains to the indexed time, we define 
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(3.5.6) 



where tj is a discrete time-point which is close to tj, so that dtj is 
smaller than At in magnitude. The local quadratic fit for a generic 
variable V(tj) is then given by [9] 


V ij (t I } " I ( d2 " d ) v iJ (t I -At) + (1-d 2 )V. j (t I ) 


+ j (d 2 +d)V (tj+At) , (3.5.7) 

where 

d = dtj/At , 

A A. 

and V^j may represent an element of either S, R, m, or fi. 

The solution for the indexed time is obtained by guessing a value 
for dtj, computing the error in satisfying (3.5.2) via 

e « ffi T (tj)<5x(t) - fi T (t I )di|) - t + tj + dtj , (3.5.8) 

and updating the value of dtj via 

dtj :« dtj - e. (3.5.9) 

Equations (3-5.8) and (3.5.9) are iteratively applied [33] until the value 
of dtj converges. The previously converged value of dtj is used as the 
starting guess in the iteration. 

The logic for propagating the gains is as follows. If dtj > At, 
then tj is incremented by At, and the gains are propagated forward by one 
time-step. If dtj < -At, then tj is decremented by At, and the gains are 
propagated backwards by one time-step. The above propagation is repeated, 
if necessary, until Idtjl < At. When tj is incremented until tj - tf, 
then the end of the maneuver is reached. Backward propagation of the 

gains, though not occur ing often, may be needed when there are 
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disturbances acting on the system, or when there is a sudden change in the 
terminal manifold during the control interval. Figure 3~2 shows the block 
diagram for the free-final-time perturbation feedback controller. 

3.6 Illustrative Examples 


The specific model considered in this section consists of a rigid 
hub with four identical elastic appendages attached symmetrically about 
the central hub, and is derived from the experimental structure of [6], 
using NASTRAN data (see Fig. 3“3). In particular, the following 
idealizations are considered: (i) single-axis maneuvers, (ii) in-plane 
motion, (iii) antisymmetric deformations, (iv) small linear flexural 
deformations, (v) only the linear time-invariant form of the equations of 
motion are considered, and (vi) the control actuator is modeled as a 
concentrated torque generating device. Figure 3 - ^ shows the first three 
antisymmetric modes, which, with the rigid body mode, defines the full- 
order model. The control system for the vehicle consists of a single 
controller in the rigid part of the structure. The structural parameters 
of the model are presented in [6]. Because of the above assumptions, only 
the antisymmetric modes are used for the example cases in this section. 
In addition, full state feedback is assumed. 

The rigid body mode and the first elastic mode are chosen for 
inclusion in the state vector for the control problem. Hence, the second 
and third elastic modes represent residual modes. Control smoothing is 
done by penalizing the first and second time derivatives of the control in 
the performance index, and augmenting the state vector with the control 
and its first time derivative. The state vector is given by 


x = [n 0 n 1 n Q u u] T , (3.6.1) 

where n 0 and n-j are the amplitudes of the rigid mode and first elastic 
mode, and u is the control torque. 
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Nominal Controls 


u N (t) 



i j 

Perturbed target 


Dashed lines indicate variables which 
depend on dt^, the change in final time. 


Figure 3-2. Block Diagram for the Free-Final-Time Perturbation 
Feedback Controller 
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As a result of many numerical simulations, it has been found that 
without modifications, the optimal perturbation feedback control, as 
presented in this chapter, performs poorly, especially towards the final 
time when the gains are large and vary quickly. One source of difficulty 
is due to the singularity of Q. Since Q is singular, all of the feedback 
gains become infinite at the final time. Because of numerical 
inaccuracies, this results in randomly large gains near the final time. A 
remedy is to use Q = -e-,1, where e 1 is a small positive number. The 
negative sign for Q is due to the fact that if Q were to be a function of 
time, it can be shown that Q < 0 for t < t f . With this modification, the 
gains become large in a well-behaved manner near the final time. 

Another difficulty manifests itself in the numerical instability of 
the final time estimation. Since the correction variable, e, of (3.5.8) 
is computed as the difference of potentially large numbers (relative to 
e), the calculation is easily numerically unstable when the values of 
m(tj) and n(tj) are large or corrupted by numerical inaccuracy. This 
results in values of e alternating in sign and increasing in magnitude at 
each successive iteration, leading to an unstable algorithm. A remedy for 
this problem is to increase the magnitude of a in the expression following 
(3.4.2) , using 



where e 2 is a small positive number. It is found that on choosing 
e 2 * 0.02, the final time estimates become much better behaved, and the 
errors in satisfaction of the terminal constraints are reduced by about 
two or three orders of magnitude. If we use e 2 = 0.5, the terminal errors 
are reduced by almost four orders of magnitude. However, with e 2 * 0*5 
the free-final-time controller behaves as if it were a fixed-final-time 
controller. From these observations, it is clear that the performance of 
the controller is extremely sensitive to the value of a. It is not 
obvious, however, that a should be modified rather than any other 
variable. Nor is it obvious what the mechanism is for the improvement of 
system performance, since the modification of a affects all of the gain 
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variables. However, it is interesting to note that the modified 
variables, Q and a, are diagonal blocks of the matrix coefficient of 
(3.4.2). The remaining diagonal block, S = Sf, also greatly affects 
system performance, as discussed in [10]. For the results of this 
section, values chosen for e., and Z 2 are 10 10 and 0.02, respectively. 


Cases 1 and 2 represent retargeting maneuvers, where the final hub 
orientation, angular rate, angular acceleration, and third time derivative 
of the hub angle are required to match a moving target whose motion is 
presumably known. The nominal target motion is a linear fly-by (see Fig. 
3 _ 5), where the target travels in a straight line at constant velocity. 
The structure is assumed to rotate about the z axis, with the appendages 
moving in the x-y plane. The components of the i|» vector are given in the 
following equations: 

0(t f ) - e T (t f ) = 0, 0 T (t f ) = tan” 1 (^) , 


0(t f ) - e T (t f ) 


0, 0 T (t f ) = 


vx 


2 2 ’ 
x + y 


0(t f ) - 0 T (t f ) 


0, © T (t f ) 


2 

-2v xy 


f 2 A 2.2 ’ 

(x + y ) 


and 


•** »«« 

0(t f ) - 0 T (t f ) 


0, © T (t f ) 


8v 3 xy 2 
(x 2 + y 2 ) 3 


3 

2v°x 


,2 2,2 

(x + y ) 


(3.6.4) 


where x = x Q and y - y Q + vt^. The perturbed target motion is also a 
constant-velocity linear fly-by, but with a different starting location 
and different velocity. The target parameters and weighting matrices are 
shown in Table 3 - 1 • The weight on the elapsed time is adjusted so that 
the nominal final time is 5.0 s, for convenience. The terminal state 
weighting is computed using a modified version of the algorithm in [10]. 
Because of the large terminal weights on n-| and , the final values of 
the remaining state variables are more or less constrained via the 
terminal constraint conditions of (3*6.3). Therefore, the elements in the 
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Figure 3-5 • Linear Fly-By for Cases 1 and 2 
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Table 3~1. Maneuver Specifications for Cases 1 and 2 


Weights 



‘ 2.63(3) 




Symmetric 


3.95(3) 

3.09(6) 





7.90 

-7. 3K2) 

2.63(3) 



s r * 






I 

1.18(1) 

4.12(4) 

3.95(3) 

2.68(5) 



6. 24(~3) 

-3.38(1) 

4.16 

1.98(2) 

1.00(5) 


_ 4. 73 (-6) 

-3.52(-2) 

4.23(-3) 

2.01 (-1 ) 

1.55(2) 3.23(3) _ 


W xx - Block diagonal [W n , W 22 , 1 . 00 ( ~5 ) , 1.00(-5)] 



2.63(“5) 

3 • 95 ( —5 )~ 


2.63(-5) 

3.95(-5)' 

" 

3.95 (-5) 

6.91 (-5)_ 

w 22 " 

_3 . 95 ( -5 ) 

6.91 ( -5 )_ 


w uu - [10] 


Nominal conditions: 

0(t o ) - 0(t o ) - n 1 (t Q ) - n-, ( t Q ) - u(t Q ) - u(t 0 ) - 0 

x 0 - 2.7(7), y 0 - -4.0(4), v - 8.0(3) 

Perturbed Conditions: 

9(t 0 ) - 2(-5)rad , 0(t o ) - -2(-5) rad/s 
ni(t 0 ) - ^(tQ) - u(t Q ) - u(t Q ) - 0 
x 0 - 2.5(7), y 0 - -4.1(4), v - 7.8(3) 

a(n) indicates a x I0 n 
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rows and columns of S f corresponding to n 0 » n 0 , u, and u can be set to 
very small numbers. However, because of numerical considerations, these 
elements are set to larger numbers to decrease the condition number of the 
Z(t) matrix. The perturbed trajectory information for Case 2 is 
introduced into the feedback law at t = 2.5 s. This simulates the case 
when target information is updated during a maneuver. 

The time history plots for Cases 1 and 2 are shown in Figures 3-6 
through 3~8. For Case 1, note that all the terminal constraints are 
satisfied for both the nominal and perturbed trajectories. (The target 
angular acceleration, 0 T , and third angular rate, 0* T , are virtually zero 
at the final time for both the nominal and perturbed targets.) For Case 
2, the hub angular acceleration shows a small terminal error. However, 
all other states reach their desired terminal values, including the 
control and its first derivative. The second derivative of the control 
torque shows a very small spike near the final time in Figure 3-8. Such 
spikes are typical when the final gains are large, and may be removed by 
appropriately adjusting S^, Q, or a of (3.4.2). 

3.7 Extensions 


The closed-form solutions given in this paper are only applicable 
for control problems where the plant dynamics is linear time-invariant. 
For nonlinear systems, closed-form solutions are much more difficult to 
obtain. Nevertheless, one may use the following method with the closed- 
form solutions of this section to approximate the solution of the feedback 
gains for a nonlinear system. First, one obtains the nominal state and 
control time-histories, using numerical techniques such as shooting 
methods [23,30] and boundary-value continuation (see Sections 4.3.1 and 
4.3.5). Second, the state differential equations are linearized about the 
nominal solution at discrete points in time. Third, piecewise linear 
time-invariant intervals are defined about these discrete points in time. 
Fourth, assuming that the feedback gains are continuous at the boundaries 
of these intervals, one can transfer the terminal conditions for S(t), 

A A A 

R(t), m(t), and n(t) into initial conditions by sequentially computing 

the boundary conditions at each of the interval boundaries, using the 
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Figure 3 - 7 . Case 2. Retargeting Maneuver with Changing 
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Figure 3~8. Control Trajectories for Cases 1 and 2 
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closed-form solutions, starting from the interval nearest the final time 
and going backwards in time. Finally, one can use the initial conditions 
for S(t), R(t), m(t), and n(t), with closed-form solutions and piece- 
wise linear time- invariant system differential equations to apply the 
feedback gains to the perturbed states and perturbed terminal constraints. 
The length of the linear time-invariant intervals may need to be adjusted 
depending on the degree of nonlinearity present at a given time along the 
nominal trajectory. 
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SECTION 4 


NONLINEAR THREE-AXIS MANEUVERS FOR FLEXIBLE SPACECRAFT 
WITH CONTROL SMOOTHING 


4.1 Introduction 


This section presents formulations for general nonlinear three-axis 
slewing maneuvers for flexible spacecraft. The approach used here is to 
find the optimal solution for the rigid body model, and then to apply this 
open-loop rigid body optimal control to the fully flexible spacecraft with 
a perturbation feedback controller. The perturbation feedback controller 
controls several flexible modes in addition to the rigid body modes, and 
the feedback gains are computed using the flexible plant linearized about 
the rigid body nominal solution at several points along the maneuver. An 
extended Kalman filter is implemented to estimate the plant states. 
Example maneuvers are shown using the model of a generic space vehicle. 

Section 4.2 presents a discussion of model development and 
simulation issues. Section 4.3 presents the solution to the nonlinear 
rigid body problem. The flexible body perturbation formulation is 
developed in Section 4.4, and the extended Kalman filter is discussed in 
Section 4.5. 

4.2 Model Development 


The spacecraft model used for the example maneuvers of this section 
is based on a satellite model similar to the N-ROSS satellite, which 
consists of a more or less rigid bus and several flexible appendages 
(Figure 4-1). For this study, the spacecraft bus is assumed to be rigid, 
and only two of the appendages, namely the radiometer and the solar array, 
are assumed to be flexible. The frequencies and mode shapes of the 
flexible appendages are in the form of NASTRAN output data. 




The original spacecraft design has one rigid body and six flexible 

appendages. Since only two of the appendages are considered flexible for 

this study, the remaining appendages and the rigid body are lumped 

together to form one rigid body. The total system inertia matrix about 

2 

the system center of mass is given in units of slug-ft by 


'3888 

-468.7 

590.7' 


-468.7 

4242 

570.2 

(4.2.1) 

_ 590.7 

570.2 

2105 _ 



The flexible appendages are each assumed to have five elastic degrees of 
freedom, and their frequencies are listed in Table 4-1 . Every mode is 
assumed to have 0.1$ damping. 

4.2.1 Multibody Dynamics Simulation 

The numerical simulation for the example maneuvers was carried out 
using a program called DISCOS (Dynamic Interaction Simulation of Controls 
and Structures) [4], which is a well-known package of software developed 
for the National Aeronautics and Space Administration (NASA) and 
distributed by Computer Software Management and Information Center 
(COSMIC). In DISCOS, a complex structure may be modeled as several rigid 
or flexible structures connected together at specific points, called 
hinges. The equations of motion for each body may then be written in the 
same general form for a single body, with the coupling between bodies 
provided by Lagrange multipliers which maintain the desired interface 
constraints. 

4.2.2 Recent Issues in Multibody Dynamics Simulation 

Doubts have recently been cast on the validity of multibody 
computer programs such as DISCOS [4], NBOD [12], ALLFLEX [15,16] and 
TREETOPS [32]. Since DISCOS was chosen to simulate the three-dimensional 
nonlinear slews of this section, an investigation was carried out to 
determine the validity of such claims [19]. 
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One of the claims was that current multibody computer programs do 
not include rotational inertia terms for the individual elements of a 
finite element model. However, as shown in Reference 9, by computing the 
inertia matrix for a single body from first principles, and comparing the 
result with the documentation for DISCOS [4], it is found that DISCOS does 
include rotational inertia terms for the individual elements of a finite- 
element model. 

Another issue of concern is the absence of a gyroscopic stiffening/ 
arc length correction term in current multibody computer programs. The 
arc length correction involves the difference in distance to a point on a 
long slender beam when measured along the deformed beam and when measured 
as the projection onto the beam's undeformed position. Inclusion of this 
correction term leads to a stiffening term in the equations of motion, 
which increases with higher angular velocities, hence the term gyroscopic 
stiffening. Since this effect is applicable only to long slender rods, 
and is important only for high angular velocities, it is often not 
included in general multibody computer programs. For the example cases of 
this study, the angular velocites are small, hence the neglected 
gyroscopic stiffening/arc length correction terms are not important. 

4.3 Optimal Nonlinear Three-Dimensional Maneuvers with Control 

Smoothing for Rigid Structures 

This subsection deals with the solution for the open-loop nominal 
control profile which is based on a rigid body model. The solution to the 
nonlinear optimal control problem is obtained by first solving the problem 
of a single-axis maneuver with a diagonal inertia matrix, and then using a 
continuation method to introduce the three-axis boundary conditions and 
off-diagonal elements of the inertia matrix. 
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4.3*1 Continuation Method 


The continuation method [29], also known as homotopy chain method, 
is a process by which a continuation parameter, a, is imbedded into the 
equations of a problem so that when a is set to zero, the modified problem 
becomes easy to solve, and when a is set to one, the modified problem 
reverts to the original hard-to-solve problem. Typically, the 
continuation parameter is used to multiply terms which make the problem 
difficult to solve. 

There are many ways of sweeping the value of the continuation 
parameter from zero to one. One way is by numerical integration which 
requires the calculation of the derivative of the solution with respect to 
the continuation. For some problems, is is difficult to compute this 
derivative. Instead, a prediction-correction type of integration may be 
performed, where previously converged intermediate solutions are used to 
estimate the required derivative by means of finite differences. A more 
crude but simple approach is to slowly increment the value of a, and 
perform an iterative correction at each increment. This approach is less 
efficient, but involves the least amount of programming. 

4.3.2 Equations of Motion 


For the rigid body control problem, let us select as state 
variables Euler parameters, f5, body angular velocities, to, pseudo- 
controls, u Q , and pseudo-control rates, u^ . The pseudo-control vector is 
defined as 


u 

o 




(4.3.1) 
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where 


'*'xy "'"xz 

I -I , 

yy yz 

I I 

zy zz -I 

and u x , Uy, and u z are torques about the x, y, and z axes, respectively. 

The torques are assumed to be applied by concentrated torque generating 
devices acting on the rigid spacecraft bus. The inclusion of the pseudo- 
controls and pseudo-control rates is for control smoothing, as described 
in Reference 9. Pseudo-controls are used instead of the actual applied 
torques because the use of applied torques results in large values for the 
angular velocity costates when the moments of inertia are large. By using 
pseudo-controls, the problem is normalized so that the values of the 
states and costates are close to the same order of magnitude. The 
equations of motion can be shown to be 


[I] - 


xx 


yx 


-i 


zx 


0 = n(w)e , (4 x i) (4.3.2) 

<*> “ u 0 + G (<*)) * ( 3 x 1 ) (4.3.3) 

u q = u 1 , (3 x 1) (4.3.4) 

and 

^ - u 2 , (3 x D (4.3.5) 

where 
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In the 

above equations, 

(4. 

3.2) is the kinematic equation 

relating the 

Euler 

parameter rates 

and 

the body angular velocities, 

and (4.3.3) 


represents Euler's equation in terms of pseudo-controls. 

4.3.3 Optimal Control Problem and Necessary Conditions 

For rigid body nonlinear three-dimensional slews, let us define the 
optimal control problem as the minimization of a finite-time quadratic 
performance index 


J = 


f^ f Cx T (t) ul(t)] W - 

' x ( t ) ' 

► 

0 

s U 2 (t) , 


h dt , 


(4.3.6) 


where 


T t T T-.T 
l8 a) u q u.J , 
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and 


W 


0 0 0 0 


0 1 


0 Q 0 0 0 

0 0 10 I/ug 


0 0 0 0 0 
0 0 I/cog 0 l/0)g 


subject to the state dynamics equations, (4.3.2) through (4.3.5), with 
specified initial and final states. The symbol I represents the (3 * 3) 
identity matrix. In the performance index of (4.3.6), the weight matrix, 
W, does not include penalties on the Euler parameters, since the angular 
displacements may be large. A weighting matrix is placed on the angular 
velocity terms so that the angular velocities may be kept small. The 
penalty terms on the pseudo-controls and pseudo-control rates are the 
time-domain equivalent of frequency-domain penalties on the pseudo- 
control, where the frequency range above oig is penalized (see Ref. 9). 


The Hamiltonian for the performance index of (4.3.6) and the state 
dynamics of (4.3.2) through (4.3.5) may be written as 


„ 1 r T„ T T . 4 T . 2n 

H = —[a) Qu) + u u + u„u_/a) D + 2u u_/a) n J 
2 OO 22B O 2 B 

+ Y^ft(a))B + ^( u 0 + G(w)) + ^0 U i + y^ u 2 > (4.3-7) 

where Y (4x1), \ (3*1), p 0 (3*1), and p-j (3*1) are costate or 
adjoint variables for 6, to, u Q , and u-| , respectively. The necessary 
conditions can then be derived from the Hamiltonian as 
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O(u)0 , 


(4.3.8) 


(i) = + G(u) , 


(4.3.9) 


r 3H iT 

L-stt-J - u . 


(4.3.10) 


r 3HiT 

hnH = u i 


U 2 ’ 


(4.3.11 ) 


[|f] T = Y = Q(»)Y , 


(4.3.12) 


[M] T - x - - o. - [£<tW - [{g] T * . 


(4.3.13) 


3H iT 


r— i 

L 3u 1 
o 


u o - U 2 /U) B - X , 


3H iT 




(4.3.14) 


(4.3.15) 


3H iT 




U 2 /u) B + U o /w B + y 1 ’ 


(4.3.16) 


where the following notation has been used 


[— ] 

L 3w J ij 


3w. * 
J 


for general vectors (or scalars) v and w. In the equations above, the 
initial and final values for $, <u, u Q and u.j are specified. However, no 
boundary conditions are known for Y, X, p Q , and . 
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4.3.4 Starting Guess for the Continuation Method 


The necessary conditions for the rigid body slewing problem shown 
in (4.3.8) through (4.3.16) represent a set of difficult nonlinear 
differential equations with split boundary conditions. To solve the 
differential equations, a continuation method is used, as discussed in 
Section 4.3.1. The first step is to find a starting guess which is easy 
to solve. For this, we choose a single-axis maneuver about a principle 
axis. To further simplify the calculations, we assume initially that the 
inertia matrix of the spacecraft is diagonal, with the non-zero off- 
diagonal terms introduced during the continuation process. 

A reasonable choice of axis for the starting guess solution is to 
use the axis with the highest peak angular momentum if the single-axis 
maneuver were to be accomplished via bang-bang control. This is also the 
axis about which the largest bang-bang torque would be applied. Denoting 
this axis by k, where k = 1, 2, or 3 corresponds to the x, y, or z axis, 
one may write a modified subset of the necessary conditions of (4.3.8) 
through (4.3.16) for rotation about axis k: 

X(t) = [K] x ( t ) , (4.3.17) 


where 


x(t) - 


*k (t) - W 


“k (t) 


[u o (t)] k 


Cu 1 (t):1 k 


C “o <t)3 k 




5(t) 
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r o i 


0 0 0 0 0 


0 


0 0 1 


0 0 0 0 0 


0 0 0 1 


0 0 0 0 


[K] 


0 

I 

k 

0 

0 

0 


. 2 
0 -0) B 0 

0 0 0 


0 -Wg 0 

0 4 -1 


0 

0 


0 0 0 -1 
0 0 0 0 

0 0 0 0 


0 

0 

0 


0 

0 

0 


0 

J_ 

2 

0 


5(t) = W - j (t-t 0 ) c, 

C = Y k (t 0 )sec(* k (t 0 )/2) , 


and <j> k (t) is the angular displacement about axis k. In deriving (4.3.17), 
we have imposed the orthogonality constraint: 


4 

l f3.Y. - 0 (4.3.18) 

, l l 


to ensure uniqueness of the Euler parameter costates. 

Since (4.3.17) represents a linear time-invariant system, the final 
values of x(t) can be related to the initial values by an equation of the 
form 


x(t f ) 


K( vv . 

e x(t G ) 


(4.3.19) 


Observing that the initial and final values for <|> k , u) k , (u Q ) k , and (u-|) k 
are known, one can perform the partitioned matrix multiplication in 
(4.3.19) for the upper partition of x(tf), and solve for the unknown 
initial values of (p 0 ) k , VM and the unknown constant C, via 


- 53 - 


' ^oW ' 


♦k (t f ) "*k (t o ) ] 


0 

( V k ( V 


Oh 

4-> 

V / 

3 

k - V 



► — 1 

'“oW 

<“oW 

c 


. <V k (t f> . 


4-> 

' 

=r 


(4.3.20) 

where 

K(t f -t 0 ) 


The initial values of B 0 , B k , Y , Y k , and X k are then given in 


terms of elements of x(t Q ) by 

B 0 (t 0 ) = cos t Q )/2) , (4.3.21) 

3 k (t 0 ) = sin (4> k (t Q )/2) , (4.3.22) 

Y o (t o ) = - C sin U k (t 0 )/2) , (4.3.23) 

\(t Q ) = C cos (<fr k (t 0 )/2) , (4.3.24) 

and 

W = 5(t 0 ) . (4.3.25) 


The initial conditions of (4.3.21) through (4.3.25), together with 
w k (t 0 ), Cu Q (t 0 )] k , Cu 1 (t 0 )] k , [y 0 (t 0 )] k , and [u 1 (t Q )] k of x(t 0 ) comprise a 
complete set of initial conditions for a single-axis maneuver with a 
diagonal inertia matrix. 


'11 


L 21 


12 


22 J 
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4.3.5 Continuation for Inertia Matrix and Boundary Conditions 


Given the single-axis rotation with diagonal inertia matrix 
starting guess of Section 4.3.4, the three-axis optimal maneuver with 
fully populated inertia matrix may be obtained through a continuation 
process. The continuation approach used is one where the continuation 
parameter, a, is incremented at discrete steps, with convergence at each 
step achieved via a Newton-Raphson iteration. 


Since there are two different quantities introduced during the 
continuation process — off-diagonal elements of the inertia matrix and 
three-axis boundary conditions — the continuation process may be performed 
separately or combined together in one process. When one combined 
continuation process is used, it may be advantageous to retain the ability 
to use separate continuation parameter increments for the two quantities 
when handling extremely difficult problems. 


For the inertia matrix continuation, the inertia matrix of (4.3.1) 
is replaced by 




I 

-a. I 
1 xy 

-a. I 
1 xz 





XX 



[K ct 1 ) ] 

- 

-a 1 I yx 

I 

-a. I 
1 yz 

• 

(4.3.26) 


yy 




- a 1 *zx 

-a. I 
1 zy 

I 





zz ^ 



Setting a 1 = 0 

produces the 

diagonal 

inertia matrix used 

in Section 4.3.4, 


and setting a-| = 1 produces the original fully populated inertia matrix. 


For the 
terminal Euler 


boundary condition continuation, let us define the modified 
angles as 


<iyt f ,a 2 ) 


<t ’j (t f ) ’ 

' a 2 <,) j (t f ) • 


j = k 
j * k 


(4.3.27) 


where k represents the axis used for the starting guess of Section 4.3.4, 
and <j>j(t f ), (j-1,2,3), are the desired final Euler angles of the three- 
axis maneuver. For each value of the continuation parameter, 02 » the 
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I 


modified final Euler parameters are computed from the values of \l> j . The 
modified final conditions for the angular velocities, pseudo-controls, and 
pseudo-control rates are similarly defined as 


uk (t f ,a 2 ) 


“j'V 

a 2 UJ j ( tf ) 


j = k , 
j * k , 


(4.3.28) 


(u o ) J <t f ,a 2 ) 


“2 (u o ) j (t f > 


j = k , 
j * k , 


and 


(u 1 )j(t f ,« 2 ) 


(u^ (t f ) 

« 2 (u 1 ) j (t f ) 


j = k , 
j * k . 


(4.3.29) 


(4.3.30) 


The modified initial conditions are defined in the same manner as for the 
modified final conditions. 


After each increase of the continuation parameters, the previously 
converged values of the initial costates no longer generate trajectories 
which satisfy the final boundary conditions. As a result, an iterative 
correction scheme is needed to correct the initial costates based on the 
error in satisfying the final boundary conditions. For this purpose, a 
Newton-Raphson first-order correction scheme is used. This is 
accomplished by Taylor expanding the terminal values of the states as 
functions of the initial costates. As a result, the partial of the final 
states with respect to the initial costates must be computed. To obtain 
quicker convergence, one may use extrapolated values of the initial 
costates based on previously converged values and back a values [8,29,34]. 

For each iteration, the modified state-costate vector is integrated 
from t Q to tf, and the error in satisfaction of the modified final 
conditions is computed. To obtain the partial of the final states with 
respect to the initial costates, one must integrate the partials of the 
state-costate vector with respect to the initial costates, along with the 
integration of the state-costate vector. That is, the matrices 
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hill 

3A(t Q ) 


and 


3A(t) 

3A(t ) 
o 


are integrated from t Q to t^., where 


A 


T t T T,T 
U A m q y 1 ] 


(4.3.31) 


and x is defined following (4.3.6). The initial conditions and 
differential equations for integrating these partials are presented in 
Reference 9. The corrected costates for the next iteration are obtained 
as follows: 


A I+1 (t o ,a) 


W a) " 


* 1 [x 


j(tj,,oc) — 


j( t t a ) ] 


(4.3.32) 


where 
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» 
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[o b 1 e 2 b 3 
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to u u. J . . f 

o 1 d tj, f a 


X I (t f ,a) = 

[0 B 1 e 2 B 3 
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0 ) u„ u, J T . , 
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In (4.3.2), the subscript <3 for x<}(t f ,a) indicates the desired modified 
final values, and the subscript I for xj(tpa) indicates the integrated 
final values for the current iterate. The variables x^Ctf* 01 )* Xi(tf» a )> 
and i represent modified forms of x<j(tf>a), Xj(tf.“)» and 
3x(t f ,a)/3A(t 0 ,a) , where the modification has been performed by replacing 
the first row of the Taylor expansion for x(tf»“) by the orthogonality 
constraint [ 36 ] 

B T (t o ,a)Y d (t o ,a) = 0 , (4.3.33) 

where Y d (t Q ,a) is the desired initial Euler parameter costate. 

The entire continuation process is summarized in algorithmic form 
as follows. (The single-axis diagonal inertia matrix starting guess of 
Section 4.3.4 is assumed to have been computed.) 


Step 1. If « 1 and a 2 ■ 1 , stop, (end of continuation). Otherwise, 
increment a 1 and a 2 , and compute [l(a-|)], x(t 0 »<* 2 ), and x(tf><* 2 ). 

Step 2. Integrate state-costate differential equations ((4.3.8) through 
(4.3.15)), and state-costate partials with respect to initial 
costates (Reference 9). 


Step 3. 

it 

Step 4. 


Compute error in satisfaction of modified final 
conditions. If small, go to Step 1. 

Compute new initial costates ((4.3.32)). Go to Step 2. 


boundary 


4.3.6 Numerical Results 

A 60 second rest-to-rest maneuver with angular displacements of 1 
radian about each axis (using a 1-2-3 Euler sequence) is simulated. The 
weighting matrix for the angular velocity is arbitrarily chosen as 
Q = 10~^I, where I is the (3 * 3) identity matrix. In choosing the value 
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1 



for the break frequency, mg, it is found that it is best to choose u)g so 
that it corresponds to the frequency of the maneuver; that is, 


2tt 


(t f" t o ) 


(4.3.34) 


For the above value of tog, the resulting maneuver has pseudo-controls with 
smooth profiles (Figure 4-2). For higher values of io B , the pseudo-control 
profiles of the resulting maneuver have more undulations (see Figure 4-3). 
This reflects the higher frequency content of the controls, directly 
resulting from the higher value of tog. For lower values of w B , the 
resulting trajectories are similar to the case where ojg is chosen 
according to (4.3.34). However, the number of Newton-Raphson iterations 
required for convergence is increased slightly, indicating that the 
partial derivative matrix of (4.3*32) may have become numerically stiffer. 
To illustrate the effect of the choice of oig on the frequency content of 
the resulting control, Figure 4-4 shows the frequency spectra of the 
pseudo-control for the single-axis starting guesses with iog corresponding 
to Figures 4-2, and 4-3. Because the penalty function of (4.3.6) 
penalizes the frequency content of the pseudo-controls for values of 
frequency above oj b , one sees a sharp roll-off near io B in Figure 4-4. 


4.4 Perturbation Feedback for Controlling the Flexible Body Response 


4.4.1 Plant Linearization and Gain Calculation 


This section presents a perturbation feedback scheme for control- 
ling the elastic deformations of a flexible body when subjected to the 
nominal rigid body torque profile of Section 4.3. The flexible plant 
dynamics is linearized about the rigid body nominal solution at several 
points in time. Steady-state feedback gains are computed based on these 
linearized plants and an infinite-time performance index with control-rate 
penalty. 

The state perturbations used for the flexible plant model for the 
perturbation feedback are 
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(4.4.1) 


- U» x 6» y «» z il cj 5* 4 S *X s *y s *z :T ■ 

where 6a) v , 5w„, and 6w„ are the perturbed body angular velocities of the 
x y z 

spacecraft bus, 5 S and 5 R are the modal amplitudes of the solar array and 
the radiometer, respectively, and 6<|> x , 6<f>y, and 6<f> z are angular 
displacements of the spacecraft bus. The modal amplitudes, £$ and 5 r* 
represent a reduced order elastic model. The numerical simulation 
includes several additional elastic modes to represent the residual mode 
responses. 

From DISCOS, the linearized system dynamics and control influence 
matrices are obtained numerically through a quadratic finite difference 
approximation. The linearized differential equation for the state 
perturbations is then 

<$x(t) = A (l) <5x(t) + B (l) 6u Q (t) , t = , (4.4.2) 

where t^ is the instant in time at which the plant is linearized, and 

B^ are the linearized state dynamics and control influence matrices at 
t = t^, and u Q (t) is the pseudo-control. 

For control-smoothing, the above differential equations are augmented by 
the differential equations for the perturbed pseudo-control : 

6u = 6u. , (4.4.3) 

o 1 

and 

6u 1 = 6 u 2 . (4.4.4) 

The performance index used for computing each set of steady-state 
gains is 
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6x(t) 


J = j [6x T (t) 6u*(t) Suj(t) fiu^t)] W 


6u (t) 
o 

fillet) 

6u 2 (t) 


y dt , (4.1*. 5) 


where 


Q 0 0 0 

0 10 I/u»g 

0 0 0 0 

2 4 

0 I/u>g 0 I/u 


The penalties on the pseudo-controls, and their rates are in the form used 
in Gupta's frequency-shaped control smoothing where the break frequency, 
u) B , may be different from the one used for the rigid body nominal. For 
each linearized plant, steady-state gains are computed based on the 
performance index of (4.4.5), and the dynamics equations of (4.4.2) 
through (4.4.4). Since the performance index is not rigorously minimized 
for the nonlinear plant, this feedback approach is suboptimal with respect 
to the performance index of (4.4.5). 

For a given plant, described by and the optimal steady- 
state feedback solution for minimizing the performance index of (4.4.5) is 
given by 


6u,(t) = - R ’[bV^ + N T ] 6 x C t ) , (4.4.6) 

d. SS 

where 

R = I/u>g , N T - [0 I/Ug 0] , 
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r A (i) 

0 
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I 0 

0 0 
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(i) 0 

0 I 

0 0 


w 


N 

R 


§ T = [0 0 I] , 


0 = P^A + aV^ - [P^B + N]R -1 [N T + bV^] + Q , 


ss 


ss 


ss 


and 


6x(t) = [Sx T (t) 6u^(t) 6u|(t)] 


During the perturbation feedback, the feedback gains are linearly 
interpolated between the points in time at which the gains are computed. 


4.4.2 Numerical Results 


Example cases are generated with the assumption of perfect state 
estimation. (Section 4.5 shows example cases where the Kalman filter is 
used for state estimation.) The 60 second rest-to-rest maneuver discussed 
in Section 4.3.6 is used for the nominal trajectory. The flexible plant 
is linearized about the rigid body nominal solution at 12 second 
intervals. Several off-nominal cases are studied. For all cases, the two 
lowest solar array modes and the two lowest radiometer modes are chosen 
for inclusion in the feedback formulation. The other higher frequency 
modes represent residual modes. All modes are assumed to have 0.156 
damping. The break frequency for the perturbation controller, ojg, is 
chosen to be half the frequency of the highest controlled mode, so as to 
minimize the excitation of the residual modes. Figure 4-5 shows the 
frequency spectra of the pseudo-control corrections when the perturbation 
feedback controller is subjected to initial conditions. 
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Figure *4-5. Typical Frequency Spectra of the 
Pseudo-Control Corrections 
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Case 1 (Figure 4-6) is the 'nominal' flexible body case, with 
perfect plant knowledge and nominal initial conditions. The controlled 
modal amplitudes and residual modal amplitudes for the solar array 
(denoted as S/A) and the radiometer are plotted separately. Note that all 
the modal amplitudes are very small by the end of the maneuver. The angle 
errors take a slightly longer time to settle. The pseudo-control 
corrections (denoted by Del u), take about 20 seconds beyond the maneuver 
time to damp out. For the case where the moments of inertia of the rigid 
bus are altered by 10$, the modal amplitude profiles are almost identical 
to those of Case 1 , while the angle error histories and pseudo-control 
corrections are altered in amplitude. 

Case 2 (Figures 4-7 and 4-8) is the same as Case 1 , except that 
initial angular errors are specified for Case 2. The error in the initial 
angle is chosen to be 5$ of the total angular displacement about each 
Euler axis. The sign of each of the errors is arbitrarily assigned. Note 
that the initial angular errors are an order of magnitude higher than the 
peak angular errors shown in Case 1 . After about 20 seconds, the angular 
error and modal amplitude time histories approach the general shapes of 
the corresponding plots for Case 1 . Since the peak values for the pseudo- 
control corrections are more than an order of magnitude higher than those 
for Case 1 , the peak modal amplitudes are also higher than in Case 1 . It 
appears that the oscillations in the pseudo-control corrections near the 
initial time may have excited the third radiometer mode (a residual mode). 
Additional example cases are presented in Reference 9. 

4 . 5 Kalman Filter for Observing the System States 

4.5.1 Gain Calculation 


This section presents a modified Kalman filter to estimate the 
system states used in the perturbation feedback scheme of Section 4.4. 
The approach presented here involves the use of linearized plant equations 
similar to those used for the perturbation feedback. 
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Figure 4-7. Case 2. Off-Nominal Initial Angles 
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Figure 4-8. Pseudo-Control Corrections for Case 2 
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Let us assume that in addition to the system dynamics matrix, 
and the control influence matrix, B^, of (4.4.2), the measurement 
influence matrix, is also linearized about the nominal trajectory at 

several points in time. For the calculation of the Kalman gains, let us 
assume that the linearized plant dynamics and measurements are subjected 
to Gaussian white noise disturbances: 


x ( t ) = A (i) x(t) + B (i) U Q (t) + w(t) , 


t = t. , (4.5.1) 


and 


y(t) = C^x(t) + v(t) , 


(4.5.2) 


where 


' 

' w(t 1 ) ' 


' w(t 2 ) ' 

T " 

' 

4 

► < 


> 

= 




/^N 

CM 

> 




Q 

0 


0 

R 


6(t 2 -t 1 ) , 


E[w(t) ] = 0 , and E[v(t)] = 0 . 


Let us assume a linear estimator of the form 


x(t) = A (l) x(t) + B (l) U Q (t) + K (l) [y(t) - C (l) x(t)] , (4.5.3) 

where is a set of constant observer gains. 

It can be shown that the gain matrix which minimizes the error 
covariance is given by 


K (i) = y( i ) [ C ( i ) ]T r 1 
ss 


(4.5.4) 


where X ^ is the steady-state error covariance matrix for the linearized 
plant: 
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0 


(it. 5. 5) 


= AX + X A T - X C T R + Q 

ss ss ss ss 


and the superscript notation has been dropped. 

During the simulation of the estimator, the variables A^\ B^, 
C^ i) , and are linearly interpolated as in Section 4.4 for the 

perturbation feedback gains. 

4.5.2 Numerical Results 


As in Section 4.4.2, the flexible plant is linearized about the 
rigid body nominal solution at 12 second intervals for the results of this 
section. The measurement variables used for the results of this section 
are the spacecraft Euler angles relative to the inertial frame; the body 
angular velocities of the spacecraft bus; the out-of-plane deformations 
and velocities of diagonally opposite corners of the solar array (points 3 
and 6); and the deformations and velocities of two points on the 
radiometer (points 9 and 10). It is assumed that raw data from sensors, 
such as accelerometers on the solar array, has already been processed to 
provide the measurements stated above. Due to the limited scope of this 
study, the sensor locations are not optimized for best performance. The 
process and measurement noise variances are chosen as small percentages of 
the peak values experienced in the 'nominal' Case 1 of Section 4.4.2. 
Case 3 (Figs. 4-9 through 4-11) shows the result of replacing the true 
state variable by the state estimate in computing the perturbation 
feedback for Case 1 of Section 4.4. Figure 4-9 shows that the Euler 
angles converge to their desired final values of 1 radian, with slight 
overshoot for the first and third Euler angles. The sensor point 
deformations shown in Figure 4-9 have very smooth profiles which converge 
to zero near the final maneuver time. Points 3 and 6 correspond to the 
two corners of the solar array, and points 9 and 10 correspond to two 
points near the center of the radiometer. The pseudo-control corrections 
of Figure 4-10 have higher peak values than the corresponding plots in 
Case 1 of Section 4.4. The angular estimate errors show the result of 
linearization at discrete points in time. One remedy is to linearize the 
plant at shorter time intervals. A better solution is to perform 
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Figure 4-11. Solar Array and Radiometer Modal Amplitudes 
and Their Estimates for Case 3 
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perturbation estimation about the nominal rigid body trajectory rather 
than estimation for the entire state. Figure 4-11 shows the amplitudes of 
the controlled modes and their estimates. One can see that the first 
solar array and radiometer modes are not estimated very well. This is due 
to observation spillover from the residual modes. In order to minimize 
the effects of observation spillover, one must choose optimal locations 
for the sensor points. Since there are 271 grid points on the radiometer 
with 5 modes, and 1000 grid points on the solar array with 4 modes (the 
original model has 15 modes), and six degrees of freedom to choose from, 
an automated procedure must be used for the selection of sensor locations 
[4,37]. However, this is beyond the scope of the current study. 
Additional examples are shown in Reference 9. 


- 76 - 



SECTION 5 


SUMMARY AND CONCLUSIONS 


This report has covered three different, yet interrelated topics. 
Section 2 has dealt with a new class of closed-form solutions for finite- 
time linear-quadratic optimal control problems. These closed-form 
solutions are used in Section 3, which presents the solution for the 
neighboring extremal path problem, as applied to spacecraft slewing 
maneuvers. Section 4 has dealt with general nonlinear slewing maneuvers 
for flexible spacecraft, for which the results of Section 3 are useful 
when the terminal conditions are slightly perturbed. A more detailed 
summary of each section follows. 

Section 2 has dealt with a new class of closed-form solutions for 
finite-time linear-quadratic optimal control problems where the plant is 
linear time-invariant. This class of closed-form solutions is based on 
Potter’s solution, which consists of a steady-state plus transient term, 
for the differential matrix Riccati equation. Five basic differential 
equations are identified for the solution of finite-time linear-quadratic 
optimal control problems. Closed-form solutions are presented for these 
five basic differential equations, and example control problems are 
presented where these solutions are used to obtain closed-form analytic 
expressions for the feedback gains, state trajectories, control 

trajectories, and residual state trajectories, with the assumptions of 
perfect plant knowledge, and perfect state estimation. 

For each example control problem, comparisons are made with closed- 
form solutions based on the Kalman-Englar method, and on the state 
transition matrix. For each case, it is found that the new class of 
closed-form solutions is more efficient than the Kalman-Englar type of 
solution based on the state transition matrix. Furthermore, it is well 
known that the Kalman-Englar solution for the Riccati matrix is 
numerically unstable when the propagation time-step is large, or when the 
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Riccati solution is not symmetrized at each time-step. Such numerical 
problems do not occur in Potter's solution for the Riccati matrix. Thus, 
it seems that the new class of solutions is numerically superior to the 
Kalman-Englar type of solutions for feedback gains and state transition 
matrix solutions for state and control trajectories. However, a rigorous 
analysis of the numerical stability and error propagation characteristics 
of the new class of closed-form solutions remains a topic for further 
research. 

The relationship between the new class of solutions and the state 
transition matrix solutions is illustrated by means of reducing subspace 
transformations for the Hamiltonian matrix. 

The closed-form solutions developed in Section 2 are applied to the 
free-final-time neighboring extremal path problem with linear terminal 
constraints, a quadratic performance index, and a linear t ime- invar iant 
plant. Closed-form solutions are presented for the perturbation feedback 
gains which cause the system to follow a neighboring extremal path when 
subjected to small perturbations in the initial conditions and terminal 
constraints. Numerical experiments indicate that slight numerical 
modifications can greatly reduce the sensitivity of the feedback gains 
near the final time. An extension is shown for using the closed-form 
solutions for problems with nonlinear plants. 

Section 4 has presented a formulation for general nonlinear slewing 
maneuvers for flexible spacecraft, whereby a rigid body nominal control 
profile is applied while a perturbation feedback controller limits the 
flexible body response and controls the plant to follow the rigid body 
nominal trajectory. The use of control smoothing in both the rigid body 
nominal solution and the perturbation feedback controller greatly reduces 
the excitation to the elastic degrees of freedom. 

Numerical results show that the break frequency used for the 
control smoothing formulation for the rigid body nominal solution should 
be linked to the maneuver time in order to produce good results. 
Numerical results for the perturbation feedback controller show that it 
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performs very well under off-nominal conditions for a 60 second maneuver. 
For further research, it is recommended that the maneuver time be 
shortened so that the break frequency for the rigid body nominal solution 
overlaps some of the structural frequencies. Such a case should prove 
challenging, since this would involve more interaction between the rigid 
modes and the elastic modes. 

A modified Kalman filter is presented for estimating the system 
states. Numerical results indicate that the approach is feasible. 
However, for further work, it is recommended that a sensor location 
optimization be performed to minimize the possibility of observation 
spillover. Furthermore, a perturbation estimation approach may be used, 
whereby one estimates the state perturbations rather than the states 
themselves. 

For maneuvers where the desired final conditions are different from 
the nominal final conditions, one must use the optimal perturbation 
feedback of Section 3. with the modifications for nonlinear plants. Such 
an approach would result in near-optimal time-varying feedback gains, 
using the same type of linearized plant as used for the results of Section 
A. 
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